Caller-held buffers and workspaces for evaluation, mass solve and contraction - #41

Caller-held buffers and workspaces for the hot paths of SimpleSplines: the one-dimensional and tensor-product spline evaluation, the KroneckerMass solve, and the tensor-product contraction and L² projection. Each new form allocates nothing; the allocating forms stay and call it. The release target is 0.4.0, a breaking release (see Breaking changes). 0.3.5, released from main, holds only the FFTA fix of #43; this branch merges it.

What changes

• evaluate!(buf, b::AbstractBSplineBasis, û, x, d = 0) and evaluate!(buf, s::Spline, x, d = 0), with a buffer of local_width(b) entries. evaluate(b, û, x, d) allocates the buffer and calls evaluate!. A wrong length raises DimensionMismatch.
• evaluate!(bufs, B::TensorProductBasis, û, x, d), with one buffer per axis; the allocating form calls it.
• mass_workspace(op::KroneckerMass, T) and mass_solve!(y, op, x, work). A workspace of the wrong element type or size raises an ArgumentError that names both values. The solve along the axes is unrolled over the tuple of axes, so an operator whose axes have different types no longer dispatches dynamically.
• contract_workspace(q::TensorProductQuadrature, T), contract!(L, q, F, work, d) and l2_projection!(û, q, F, work). The contraction is unrolled over the axes too.
• The result type of an evaluation: where the basis and the coefficients are IEEE floats (Float16, Float32, Float64) or complex numbers of them, and the point is an IEEE float, the result takes the precision of the point. Every other type is kept by promotion, as in 0.3.5. Knot differences are formed in the precision of the basis before they are narrowed, so a Float32 evaluation far from the origin stays finite and accurate.
• Tests in test/integration/workspaces.jl: == 0 bytes after a warm call through a fixed-arity barrier, beside a > 0 control, behind the coverage guard; == against the allocating form and @inferred, in Float32 and Float64, on clamped, periodic and recombined bases, 2D periodic × clamped and 3D tensor products, and ForwardDiff.Dual arguments with Dual buffers; ForwardDiff.gradient with respect to the coefficients; BigFloat and complex coefficients, also at derivative orders above zero. ForwardDiff joins the test dependencies.
• CHANGELOG.md opens ## [Unreleased] — targeting 0.4.0, with Breaking Changes and the bytes before and after.

Breaking changes (0.4.0)

• An evaluation at an IEEE float point returns the precision of the point. This changes the result wherever the point is narrower than the promotion of the basis and the coefficients: a Float64 spline at a Float32 point returns a Float32 (0.3.5: Float64), so does a Float32 basis with Float64 coefficients at a Float32 point, complex coefficients return a ComplexF32, and evaluate_all returns buffers of Float32. It holds for evaluate, evaluate_all and Spline evaluation on 1D, tensor-product, recombined and polar bases. Remedy: pass the point in the precision of promote_type(eltype(basis), eltype(û)), for example Float64(x) or Float64.(x). Dual points or coefficients still give a Dual, BigFloat gives a BigFloat, and Int and Rational points give the basis type.
• mass_solve!(y, op::KroneckerMass, x) and ldiv!(y, op, x) raise DimensionMismatch for an x whose length is not length(y). In 0.3.5 a shorter x solved on stale entries of y without an error, and a longer x raised a BoundsError. Remedy: pass an x of length(y). \ is not affected.
• The new caller-buffer forms refuse some buffers. They are new, so no 0.3.5 call changes: evaluate! refuses a buffer that does not hold the point type or is not floating-point (ArgumentError), and the workspace forms refuse a workspace of the wrong element type or size (ArgumentError). Remedy: zeros(typeof(x), local_width(b)), mass_workspace(op, eltype(y)), contract_workspace(q, eltype(L)).
• Downstream: PoissonSolvers 0.6.x admits SimpleSplines = "0.2, 0.3", so a package that uses both needs a PoissonSolvers release that admits 0.4.

Verification

• The full suite at 037f9c3 passes on Julia 1.13.3 (69 047 tests, the allocation assertions active) and on the 1.12.7 floor (69 047 tests, --check-bounds=yes), with FFTA 0.3.2 and FFTW 1.10.0. The commit after it changes only CHANGELOG.md.
• A separate critic judged the changes since the review in two rounds. The code passed: the result-type rule on 1D, tensor-product, recombined, Natural and polar bases with Float16/Float32/Float64, complex, Dual, BigFloat, Int and Rational types; ForwardDiff.gradient with respect to the coefficients; the DimensionMismatch on all four solve routes; no unbound type parameters on Julia 1.13 and 1.14; 0 bytes and @inferred on the buffer forms; every mutant of the fix caught. The last round failed on one sentence of the Breaking Changes entry, which 44ed230 corrects.
• CI on this branch: every required job passes. The Downgrade job fails as it does on main (a DocStringExtensions/Documenter resolver conflict), and it is advisory.

Open for the maintainer

• The allocating contract and three-argument l2_projection! allocate more than on 0.3.4 (29 584 / 14 912 and 29 344 / 14 880 bytes against 24 064 / 13 536 and 25 296 / 14 640, Float64 / Float32), because each builds two buffers of the largest intermediate array; the CHANGELOG states it, and no test pins it.
• _check_coefficients(û, q::TensorProductQuadrature) takes its arguments in the opposite order to _check_coefficients(basis, c) of src/spline.jl.
• The single-function forms of a recombined and of a polar basis sum parent values that are already narrowed to the point type, where evaluate_all! sums them before narrowing; the two can differ by round-off in the point's precision.
• The pole sums of src/polar.jl run in the result type; a test cannot tell that from the coefficient type (a difference of about 5 ulp).

Known issues

KNOWN_ISSUES.md adds K6 (a periodic basis gives NaN derivatives to a ForwardDiff.Dual at a and b), K7 (l2_projection! into a Float32 array fails on a Float64 periodic quadrature) and K9 (evaluate! on a TensorProductBasis takes buffers of one type only, as evaluate_all! does).

🤖 Generated with Claude Code